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Abstract 

A geometric derivation of numerical integrators for optimal control problems is 
proposed. It is based in the classical technique of generating functions adapted to the 
special features of optimal control problems. 

1 Introduction 

Optimal control has been one of the driving forces behind many of the applications 
of mathematics to engineering, robotics, economics... In fact, the Maximum Principle 
was discovered by L.S. Pontryagin in 1955 in an attempt to find a solution for a 
highly specific optimization problem related to the manoeuvres of an aircraft. One 
of its main features is the interplay among different research areas, specially control 
theory, classical mechanics and differential geometry. Historically, Optimal Control 
Theory (OCT) took place during the 1950's and its geometrization was started in the 
1960's. This geometric analysis of OCT has been introduced using many fundamental 
tools of differential geometry: Lie groups, exterior differential systems, fiber bundles, 
riemannian and subriemannian geometry among others. 

From other point of view, a geometric methodology has been recently shown to 
be very useful for simulating numerically the motion of dynamical systems. Following 
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this research line, new numerical methods have been developed, called geometric inte- 
grators; usually, these integrators, in simulations, can run for longer times with lower 
spurious effects (for instance, bad energy behavior for conservative systems) than the 
traditional (non-geometrical) ones. In particular, we are interested in extensions to 
OCT of Discrete variational integrators. These integrators have precisely their roots in 
the optimal control literature in the 1960's and 1970's (Jordan and Polack |JorPol:64j . 
Cadzow |Cadz:70j . Maeda jMESHOl EEEM]) and in 1980 ' s b Y Lee [EHEHHl EEEZ1 , 
Moser and Veselov [MosVes:91 . Although this kind of symplectic integrators have been 
considered for conservative systems LTar Nor:97a| [KaMaO r799"l IMarWes:f)lj . it has been 
recently shown how discrete variational mechanics can include forced or dissipative 
systems KMOW:00, MarWes:01 , holonomic constraints Ma rWes:01j . time-dependent 
systems (LeoMd D : 20021 1 M ar WesTfiTj . frictional contact PKMO:02] and nonholonomic 
constraints (see jCort:021 ICorMar:01l [LeMDSa:02al ILeMDSa:02bj ^. Moreover, it has 
been also discussed reduction theory |Bob Sus:99a. BobSus:99b], extension to field the- 
ories |JarNor:97bl IMa PaSh:98 and quantum mechanics |NorJar:98] . 

In this paper, we shall continue this work by extending to the discrete variational 
techniques to Optimal Control Problems and relating our results with Discrete Opti- 
mal Control Theory. Mainly, we shall give a geometrical construction of symplectic 
integrators for OCT, proving as a direct consequence the symplecticity of some discrete 
optimal control problems. As a nice consequence, an easy proof of the symplecticity 
of discrete Hamilton equations will be given. 

Since most engineering systems are time-dependent, we shall include the time vari- 
able explicitely in our control models and some geometric tools (mainly, cosymplectic 
geometry) of time-dependent mechanics will be useful 

2 Optimal control theory 

It is well known that the dynamics of a large class of engineering and economic systems 
can be expressed as a set of differential equations 

q A = T A (t,q(t),u(t)), 1 < A<n, (1) 

where t is the time, q A denote the state variables and u a , 1 < a < m, the control inputs 
to the system that must be specified. Given an initial condition of the state variables 
and given control inputs we completely know the trajectory of the state variables q(t) 
(all the functions are assumed to be at least C 2 ). 

Given an initial condition, usually qo = q(to), our aim is to find a C 2 -piecewise 
smooth curve ^{t) = (q(t),u(t)), satisfying the control equations Q and minimizing 
the functional 

J( 7 )= ! L(t,q(t),u(t))dt + S(T,q(T)), (2) 

for some fixed and given final time T G R + . The integral L(t, q(t),u(t)) dt depends 
on the time history (from to to T) of the state variables and the control inputs, and 
S(-,q(-)) is a cost function based on the final time and the final states of the system. 

In a global description, one assumes a fiber bundle structure it : IR x C — ► Q, 
where Q is the configuration manifold with local coordinates (q A ) and C is the bundle 
of controls, with coordinates (q A ,u a ), 1 < A < n, 1 < a < m. 
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The time-dependent ordinary differential equations (JJ) on Q depending on the 
parameters u can be seen as a vector field V along the projection map tt, that is, T is 
a smooth map r : WL x C — ► TQ such that the diagram 

RxC 




9 A a 

is commutative. This vector field is locally written as V = V (t, q, u) 



' dq A ' 

A neccesary condition for the solutions of such problem are provided by Pon- 
tryaguin's maximum principle. If we construct the pseudo-Hamiltonian function: 

H(t, q,p, u) = PA^ A {t, q, tt) - L(t, q, u) = pT(t, q, u) - L(t, q, u) (3) 

where pa, 1 < A < n, are now considered as Lagrange's multipliers, then a curve 
7 : [to 3^1 ~* C, 7(t) = (q(t),u(t)) is an optimal trajectory if there exist functions 
PA(t)> 1 ^ -A < n, such that they are solutions of the pseudo-Hamilton equations: 



dH 

q A (t) = —(t,q(t),p(t),u(t)) 
OPA 



(4) 



PA® = -^A^Q(t),P(t),u(t)) 



dH 

d<p 

and we have 

H(t,q(t),p(t),u(t)) = min H(t,q(t),p(t),v), t G [t ,T] (5) 



q(0) = q Q and Pa{T) = -^-(r,a;(T)) 



with transversality conditions 

OS 
dq 

Condition (jSJ) is usually replaced by 

dH 

— =0, l<a<m, (6) 

when we are looking for extremal trajectories. 

It is well known that the Pontryaguin's necessary conditions for extremality have 
a geometric interpretation in terms of presymplectic (or precosymplectic) Hamilto- 
nian systems. The total space of the system will be K x (T*Q xq C), with induced 
coordinates (t, q A ,pA, u a ) . 

Define the Pontryaguin's Hamiltonian function H : R x (T*Q xq C) — ► R as 
follows 

H(t, a q , u q ) = (a q , T(t, u q )) - L(t, u q ) 
where a q G T*Q and (t,u q ) G ir^ 1 (q). Therefore, the coordinate expression of H is 

Let loq = —<19q be the canonical symplectic form on T*Q, where 9q is the Liouville 
form, and consider the canonical projection tt\ : IR x (T*Q xq C) — > T*Q. Define 
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the 2-form Q H on E x (T*Q x Q C) by n H = tt^q + dH A dt. Then, (dt,tt H ) is a 
precosymplectic structure onRx (T*Q xq C) (see LeMaMD:96j). 
Eqs. (j3J) and (JSJ can be intrinsically written as 

i x n H = 0, i x <ft = 1 (7) 

Since (dt, is a precosymplectic structure, Eqs. (J7J) need not have a solution, in 
general. 

Applying the Dirac-Bergmann-Gotay-Nester algorithmn |Dirac:64l IGotNes:79| to 
the precosymplectic system 

(M x (T*Q x Q C),dt,Q H ,H) 

(see |ChLeMa:94j ) we obtain that Eqs. © correspond to the primary constraints for 
the precosymplectic system: 

Y du a 

Eqs. Q have algebraic solution along the first constraint submanifold Pq determined 
by the vanishing of the primary constraints. On the points of Po there is at least a 
pointwise solution of Eq. (J7J), but such solutions are not, in general, tangent to Po. 
These points must be removed leaving a subset Pi C Po (it is assumed than P\ also is a 
submanifold). Thus, we have to restrict to a submanifold P2 where the solutions of |7J) 
are tangent to Pi . Proceeding further this way, we obtain a sequence of submanifolds 

... M p^... M p 2 ^ Pl ^P oM l x (t*q XQ C ) 

If this algorithm stabilizes, i.e. there exists a positive integer k £ N such that Pfc = 
Pfc + i and dimPfc ^ 0, then we shall obtain a final submanifold Pf = P^, on which a 
vector field X exists such that 

(i x n H )\P f =0, [i x dt=l)\p s (8) 

The constraints determining Pf are known, in the control literature, as higher order 
conditions for optimality. 

If X is a solution of |HJ) then every arbitrary solution on Pf is of the form X' = X+£, 
where £ G (ker n H n ker dt) n TP f . 

Therefore, a necessary condition for optimality of the curve 7 : R -» 1 x C, 
j(t) = (t,q(t),u(t)) is the existence of a lift 7 of 7 to Pf such that 7 is an integral 
curve of a solution to Eqs. (jHJ). 

In the regular case, the final constraint manifold will be Po (that is, Po = Pf) 
and all the constraints are of the second kind following the classification of Dirac (see 
|LeM aMD:96]). In such case, (Pq,Q,t]) is a cosymplectic manifold, where Q and rj 
denote the restrictions of £Ih and dt to the submanifold Po. Denote also by uj and 9 
the restrictions of t^Iujq and n*0Q to Po. 

The cosymplecticity of (Po, rj, Q) is locally equivalent to the regularity of the matrix 

d 2 H 



du a du b 



l<a,b<m 

along p> The dynamical equations for the optimal control problem will become 

i x n = 0, i xV = 1 (9) 
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Taking coordinates (t, q A ,PA) on Pq-i then @ are equivalent to: 

dH 

q A (t) = ^(t,q(t),p(t)) 
PA(t) = -^(t,q(t),p(t)), 



(10) 



where we have substituted in (jlj) the control variables u a by its value u a = f a (t,q,p), 
applying the Implicit Function Theorem to the primary constraints (j) a = 0. This also 
implies that we have a canonical projection from Pq onto M, say 7To : Pq — > R . 
In such case, there exists a unique solution Xp of Eq. @: 

ix PQ tt = 0, i Xpo r] = l 

and its flow preserves the cosymplectic structure given by f2 and r/. That is, if we 
denote by Fh the flow of Xp then F?Cl = il, and F^ij = r]. In local coordinates, 

(2) (2) 

Fh(t ,qo,Pa) = {t + h,qi,pi). Denote by F^ the mapping F^ '{t ,q ,p ) = (qi,pi), 
and by F tl ^ '■ Pq° — > Pq 1 the mapping defined by 

Fh,t {qo,Po) = F^l to (t ,q ,p ) , 

where we write Pq = (ttq)" 1 ^), with t £ I. Obviously, Pt 2 ,ti ° F tl>to = F t2t t in their 
common domain. 

The submanifolds Pq naturally inherit a symplectic structure Ut by taking the 
restriction of a; to Pq. Similarly, denote by 0t the restriction of 6 to Pq, then u>t = —dOt- 

It is easy to deduce that, in such case, F tlt t is a symplectomorphism; that is, 
F ti,to u ti = Uto, noting that 

Q = uj + dH\p A rj 

This last remark will be interesting for constructing geometrical integrators for explic- 
itly time-dependent optimal control systems. 



3 Generating functions 

Let (Mi, Ui), i = 0, 1 be two exact symplectic manifolds (i.e. u>i is symplectic and exact, 
Ui = —d9i, i = 0, 1) and suppose that g : Mq — > M\ is a diffeomorphism. Denote by 
Graph(g) the graph of g, Graph(^) = {(xo,g(xo)) / xq G Mo} C Mq x M\. Denote by 
7Tj : Mo x Mi — »■ Mj, i = 0, 1 the canonical projections, and consider the 1-form and 
2-form on Mo x M\ defined by 

@(i,o) = "l^i - ttq9 

^(1,0) = tti^i - ^o^o = -d6(i,o) 



As it is well known ^(1,0) is a symplectic form. 

Let i g : Graph^) <^-> Mo x Mi be the inclusion map, then 

■fl n (l,o) = (^0|Graph(g))*(5*^i -^o) 
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Using this equality, it is clear that g is a symplectomorphism if and only if i*J7( 10 ) = 0, 
that is, if Graph(g) is a Lagrangian submanifold of (Mo x Mi,On )). 
Now, if g is a symplectomorphism we have 

^(1,0) = -di* g @(i,o) = 
and, therefore, at least locally, there exists a function S : Graph (g) — > IR such that 

i*e (li0) = ds (ii) 

Let (qo,Po) and (gi,pi) Darboux coordinates in Mo and Mi, respectively. Since 
Graph(g) is diffeomorphic to Mo, we can take (qo,Po) as natural coordinates in Graph(g) 
Since (qo,Po,qi,Pi) are coordinates in Mq x Mi, then, along Graph(g), we have 
qi = qi(qo,Po), Pi =Pi(qo,Pa) and 

pi dqi -p dqo = dS(q ,Po) 



3.1 Generating functions of the first kind 

Assume that in a neighborhood of some point x £ Graph(g), we can change this 
system of coordinates by new independent coordinates (qo,qi) (the local condition is 
that det (dq\ / dpo) ^ 0). In such case, the function S can be expressed locally as 
S = S(q ,p ) = Si(q ,qi). 

Definition 3.1 The function Si(qo,qi) will be called a generating function of the 

first kind of the symplectomorphism g. 



From (fTTj) we deduce that 



Po 



Pi 



dSi 

dqo 
dSi 
dqi 



(12) 



(see Fig. 1) 



Conversely, if Si(qo,qi) is a function such that det 



( a 2 Si 

\dq dqi 



^ then Si(q ,qi) 



(po,Pi) is a generating function of some canonical transformation g implicitly deter- 
mined by Eqs. (O, g(qo,Po) = (<?i,Pi) (see |Arn:78j). 

^o,-f|) 



{qo,qiY 



Fiff.l 



(^OlGraph^) 

/_ dSi _ dSi \ 

Wo,--a^,qi,-9^) 
( 7r 2)|Graph( 9 ) 
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Now suppose that M is a fiber bundle over the real line R, tt : M — > R, and 
Mt = 7r _1 (i) are the fibers, where each fiber Mf is equipped with a symplectic form 
ujf Let 3( S) i) : M< — > M s be a two-parameter family of symplectomorphisms satisfying 

9(h,t x ) 9(h,t ) = 9{t 2 ,t ) 

Next, we shall show how this composition law can be translated in terms of their 
respective generating functions. Moreover, the following results will give a geometric 
interpretation of the Discrete Euler-Lagrange equations MarWe s:01| . 

Theorem 3.2 Let 5| tjv '*°) oe a function defined by 

N-l 

st' t0 \q ,q N ) = J2 S l k+1 ' tk) (^,qk + l) 

where q k £ Mt k , 1 < k < N — 1, are stationary points of the right-hand side, that is 

= D 2 sf k ' tk - 1 \q k ^,q k )+D 1 s{ tk+1 ' tk \q k ,q k+1 ), 1 < k < N - 1. 

If sf k ' tk ~^ are generating functions of the first kind for 3(t fe ,t fc _i)^ then S t 1 N ' t °' > is a 
generating function of the first kind for g(t N ,t ) '■ Mt Q — > M tN . 

Proof: Recursively, it is suffices to give the proof for N = 2: 

St to \q ,q2) = S^ t0 \q ,x) + St U \x,q 2 ) 

where x is an stationary point of the right-hand side. 

From the definitions of generating functions for <7(t 2i j 1 ) and g^fy) 

Pldqi-podqo = dS ,(il '* o) (g , qi) 
P2 dq 2 - Pi dqi = dsf 2 M ^ (qi , q 2 ) 

and therefore 

P2dq 2 -p dq = diS^^q^+S^iquq^) 

It follows that 

= D 2 st t0 \q ,qi) + D 1 s[ t2 ' tl \q 1 ,q 2 ) 
and, obviously, for this choice of q\ then 

St(qo,qi)+SZ(qi,q 2 ) 

is a generating function of the first kind of g(t 2 ,t )- I 

Now, we are in condition to bring this procedure to the limit when the number of 
subintervals increases to infinity. Consider as its continuous counterpart a cosymplectic 
manifold (M,rj,u>), where M is still a fiber bundle over R (tt-^ : M — ► R) and ry = 
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7rjg,((fi). Denote by Mt = 7^ (t), t 6 1, Take a Hamiltonian function H : M — > IR 
and its Hamiltonian vector field Xh given by 

ix H u = and ix H f1 = 1 

Let -F^s) : M s — > be the two-parameter family of symplectomorphisms generated 
by X# (see section [^J) and consider as symplectic form on each fiber Mt the restriction 
of u to this fiber. 

We shall give a characterization of the generating functions of the first kind asso- 
ciated to Fu ;S ) for t close enough to s. For doing that, consider Darboux coordinates 

/ Q^H \ 

(t,q A ,pA) on M and assume the regularity condition det I — — - — I ^ 0. Thus, 



dp A dps 

Proposition 3.3 A generating function of the first kind for Fu s -\ is given by 



St t0) (q Q ,qi) = [ 1 (p(t)q(t)-H(t,q(t),p(t))) dl 
Jt 



where t — > (t,q(t),p(t)) is an integral curve of the Hamilton equations such that q(to) 
q and q(ti) = q\. 



Proof: We only use Hamilton equations and integration by parts 

i (qo,qi) = / i 

dqo J t \dqo 



Qgyt 

(qq,qi) = I ( T^q + PTr- - ^ir- - " Q " *^ ) dt 



and 



t 

dqi 

-Po+Pi-^— 
dqo 



dSt' t0 \ v f h {dp . . , 

-{qq,qi) = / \^—q + p^ ^-^—) dt 



dq 


dH dq 


dH dp 


dqo 


dq dq 


dp dqo 


dq " 


j dt 




dqo, 






-Po 




dq 


dH dq 


dH dp 



dqi Jt \dqi dqi dq dqi dp dq 1 

'■ ' dq ,.dq\ 

t \ dqi dqi J 

dqo _ 
z Pi-Po^— =Pi I 
dqi 



Remark 3.4 Suppose that t^i —ti = h, for all i = 0, • • • N — 1, then from Theorem 
13.21 we have 

JV-l 

S? h (qo,q N ) = J2 S l(l>°> ( lk+i) 

k=0 

where 

= D 2 S^(q k _i,q k ) + D 1 S%(q k ,q k+1 ), 1 < k < N - 1. 
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Now, if we take as new generating function an adequate approximation S^j of S± then 
= D 2 SjHq k - U q k )+D 1 S%(q k ,q k+1 ), l<k<N — l. 

are the well-known Discrete Euler-Lagrange equations (see MarWcs:01 and references 
therein). For instance, one can take 

Sd(qo,qi) = hC(aq + (1 - a)q 1: 1 ° ), a 6 [0, 1] 

or alternatively, we could have considered more accurate approximations. Here, we 
are assuming that £ : I x TQ — > R is a Lagrangian function related via Legendre 
transformation with the Hamiltonian function H (see |Arn:78j ) which is locally possible 
because of the regularity of H . 

Denote by Si(qo, qi, to, ti) = S^ 1 (qo,qi). From Proposition (j3,3[) . it is easy to 
show that: 

D 3 S 1 (q ,q 1 ,t ,t 1 )=D 3 S^' to \q ,q 1 ) = H{t ,q ,p ) 
D i S 1 (q Q ,q 1 ,t ,t 1 )=D 4 S^' to \q ,q 1 ) = -H{tx,qx, Pl ) 

(see also |MarWes:01| ). As a consequence 

D A S^ tk -^\q k - U q k ) + D 3 S^+^\q k ,q k+1 ) = (13) 

It should be noticed that if we take a new function S$ k+1,tk) as an adequate ap- 
proximation of Sf( f k+i>*k) ? then solutions {go, qi, ■ ■ ■ , qN} of equations 

D 2 S i J k ' tk - 1 \q k -i,q k ) + D 1 S$ k+utk \q k ,q k+1 ) = 0, 1 < k < N - 1. 

do not satisfy (|T3|) for arbitrary values of t k -\,t k ,t k+ \. Therefore, we may write the 
system of difference equations 

D 2 S^ tk -'\q k ^q k ) + D 1 S^ tk \q k ,q k+l ) = ^ 
D A S^ tk -i\q k ^,q k ) + D 3 S^ tk) (q k ,q k+1 ) = 0, 

which under regularity assumptions will determine a time-dependent discrete flow 

&(qk-i,qk,t k -i,t k ) = (q k ,q k+ i,t k ,t k+ i) 

with variable step size h k = t k+1 - t k (see |KaMa( )r:99l lhee:83l lhee:87l ILeoMn]P:2002l 
IMarWes:01p . 

3.2 Generating functions of the second kind 

The construction of more general generating functions will be useful in next sections. 
For instance, suppose that (go,.Pi) are independent local coordinates on Graph(g). 
Then the function S is written as S = S(qo,pi). 
We have 

pi dqi - p dq = -qi dpi + d{qi Pl ) - p dq = dS. 

If we define 

^2 (go, Pi) = qiPi ~ S{q ,Pi), 
where qi is expressed in terms of go an d P \, then we deduce that 

qidpi +Podq = dS 2 (q ,pi) 
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Definition 3.5 The function S^QOjPi) w ^ be called a generating function of the 

second kind of the symplectomorphism g. 

We have that 

( dS 2 

Po 



(<9 2 S2 \ 
— — - — ) 7^ then 
dq dpij 

5*2 is a generating function of some local symplectomorphism determined by Eqs. (|15j) 
(see |Xrni78j). 

Denote by F/ t2itl \ : — > M^ 2 the two-parametric group of canonical transfor- 
mations generated by the Hamiltonian vector field Xjj, as in the preliminaries to 
Proposition 13.31 We have the following. 

Theorem 3.6 Let a function 5^'*°) be defined by 

N-l N-l 
St ,t0 \Q0,PN) = £ Si k+1 ' tk \q k ,p k+1 ) - £ QkPk (16) 

k=0 k=l 

where q k , 1 < k < N, and p k , < k < N — 1, are stationary points of the right-hand 
side, that is 

q k = (q k -i,Pk), l<k<N, (17) 

Pk = (q k ,Pk+i), 0<k<N-l, (18) 

then S^' 10 ^ is a generating function of the second kind for Fu Njto \ : Mt — > Mt N . 

Proof: It follows as in Theorem 13.21 | 
As a consequence, we have that 

N-l 

S {tN ' to) (go, Pn) = QnPn ~ st' to \qo,PN) = £ [^k+iPk+l - S { * k+1 > tk \q k ,p k+1 ) 

k=0 

(19) 

where the unknown coordinates are given by (|17|) and (|18l) . 

Proposition 3.7 A generating function of the second kind for F^ tl to ^ is given by 

st tt,] (qo,Pi) =Piqi - / (P(t)q(t) ~ H(t,q(t),p(t))) dt 

where t — > (q(t),p(t)) is an integral curve of the Hamilton equations such that q(to) = qo 
and p{t\) = pi. 
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Proof: It is proved in a similar way to Proposition 13.31 | 

Denote by S^t, %,Pi) = ^ (QOjPx) then it is easy to show that (see, for instance 
|MaLuWa:()2j ) 

Theorem 3.8 (Hamilton-Jacobi equation for S2) // 5*2(t, qo,Pi) is a solution of 
the partial differential equation 

= H(^-(t,q ,pi),pi), S 2 (0,q ,pi) = q pi (20) 

then the mapping (qo,po) — ► (qi,Pi) defined by Eqs. M5)) is the exact flow of the 
Hamiltonian system determined by H . 

4 Optimal control of Discrete-time systems 

In this section we shall define the general solution of an optimization problem for 
discrete systems and analyze its geometric behaviour, in particular, the symplecticity. 
Suppose that the discrete state equations are given by the dynamical equation 

Qk+i = f A (k,q k ,u k ), fc = 0,l,...,JV-l, A = l,2,...,m (21) 

or, shortly, q k +\ = f(k,q k ,u k ), where qo is initially given. 
The associate performance index or objective function is: 

iV-1 

J = S(N, q(N)) + L(k, q k ,u k ) (22) 
fc=0 

where S is a function of the final time and state at the final time N, and L is time- 
varying function of the state and control input at each intermediate discrete time k. 

The optimal control problem is solved finding controls lit, k = 0, 1, . . . N — 1, that 
drive the system along a trajectory g|, k = 0, 1, . . . , N, verifying the state equations 
such that the performance index is minimized. 



4.1 Problem solution 

Let us now solve the optimal control problem for the discrete optimal problem de- 
termined by ()21j) and (|22|) using the Lagrange multiplier approach. Considering the 
state Eqs. (|21jl as constraint equations, then we have N ■ m constraints, and we as- 
sociate a Lagrange multiplier to each constraint. Next, we construct the augmented 
performance index J' by 

AT-l 

J '= 2 [Pk+l(f(k,q k ,n k ) - q k+1 ) - L(k,q k ,u k )} - S(N,q(N)) (23) 

k=0 

where p k +\ = ((pfc+i)^) are considered as Lagrange multipliers with A = 1, . . . , n and 
fe = 0,...,JV-l. 

Taking the Hamiltonian function 

H(k,q k ,p k+1 ,u k ) = p k+ if(k,q k ,u k ) - L(k,q k ,u k ) 
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we deduce that the necessary conditions for a constrained minimum are thus given by: 
dH 

q k+ i = -7^(k,q k ,p k+1 ,u k ) = f(k,q k ,u k ) (24) 

dH df .. . dL . 

Pk = -r^-(k,q k ,p k+1 ,u k ) = p k+ ± — {k,q k ,u k ) - — (k,q k ,u k ) (25) 

dH df , dL . 

= ^{k,q k ,p k+ i,u k ) =p k+1 — (k,q k ,u k ) - —{k,q k ,u k ) (26) 

where < < iV — 1, and the transversality conditions 

dS 

PN = —Q-{N,qN) and q fixed. 

Observe that the recursion for the state q k develops forward in time, but the co-state 
variable p k develops backwards in time. Therefore the required boundary conditions 
for finding a solution are the initial state go an d the final co-state p^. 

Assume that ^ _ 

det (du a dub) ^ ° 
then, locally, u* k = h{k,q k ,p k+ \). If we denote, by 

H(k,q k ,Pk+i) = H(k,q k ,Pk+i,u*k) 
then Eqs. 1)24(1 . (|23|) are rewritten as 

dH 

Qk+i = -Q^(k,q k ,p k+ i) (27) 
dH 

Pk = -r^-(k,q k ,p k +i) (28) 

with < k < Ni. 

Consider the function 

G k (q k ,q k+ i,p k+1 ) = H(k,q k ,Pk+i) -Pk+iQk+1, < k < N - 1. 



Then, for a fixed k: 

dH dH 

dG k = -—(k,q k ,p k+ i)dq k + - (k, q k ,p k+1 ) dp k+l - p k+ i dq k+1 - q k+ i dp k+ i . 

dq k dp k+1 

Along solutions of Eqs. ((HJ), (|2*5")) and ((2f)|) we have: 

^fcjSol = Pk dq k — p k+ \ dq k+ i , 

which implies 

dpk A dq k = dp k+ i A dq k+1 . (29) 

along the solution of P ^ -ipB j) . 

In the next subsection, we shall analyze the geometric meaning of Eq. (|29|) . which it 
is obviously interpreted as symplecticity of discrete optimal control problems in terms 
of a natural symplectic form. 
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4.2 Generating functions of the second kind and discrete 
optimal control problems 

From Proposition 13.31 the following function is a generating function of the second 
kind for the cosymplectic Hamiltonian system (Pq, r/, Q, H\p Q ), which determines the 
dynamics of the optimal control problem given by and JSJ): 

st t0 \qo,Pl)=Piqi- {p(t)q(t)-H lPo (t,q(t),p(t))) dt , (30) 

where t — ► (t,q(t),p(t)) is the integral curve on Pq of the vector field Xp . Here Xp 
is the unique solution of equation 

ix Po Q = dH\p i Xpo r] = l 

with (q(t ),p(t )) = (q ,p ) and {q(h) , p(h)) = (<?i,Pi). 

We now focus on the construction of a numerical integrator for the Hamiltonian sys- 
tem (Pq, r/, f2, H\p Q ) by using an approximation of the generating function. As we shall 
show, the obtained method also realize the integration steps by symplectomorphism 
transformations; then, it is a symplectic integrator. 

First take a fixed time interval h = tk+i — t/., k = 0, . . . , N — 1. 

Assume that we are working on vector spaces, and consider the following natural 
approximation : 

S%(k,q k ,p k+1 ) = Pk+iqk+i ~ hpk+i f Qh+1 , qk ) - hL(k,q k ,p k+1 ) 



h 

+hp k+1 t(k,q k ,p k+1 ) 

where, for instance, L(k,q k ,p k+1 ) = L\ Po (t +kh,q k ,p k+1 ) and f(k, q k ,p k+ i) = F| Po (t + 
kh,q k ,p k+1 ). 

If we denote by f(k,q k ,p k+ i) the function 

f(k,q k ,p k +i) = hf(k,q k ,p k+ {) + q k (31) 

then, 

S2(k,q k ,p k+ i) = p k+ if(k,q k ,p k+1 ) - L(k,q k ,p k+1 ) = H(k,q kl p k+1 ) . 
Thus, equations 



— j -(k,q k , Vk+1 ) = — k -{ 



Pk = -^Zk( k ><lk,Pk+l) = -^j:(k,q k ,p k+1 ) 

(32) 

qk+i = ^ {k,qk,Pk+i) = ^ {k,qk,Pk+i) 

OPk+l OPk+1 

are exactly ()27|l and (|28[) and the symplecticity condition (|29|) for discrete optimal 
control problems is now a trivial consequence of the generating function construction. 

Remark 4.1 It is also possible to construct symplectic numerical methods of higher 
order; for instance, considering better approximations of the Hamilton Jacobi equation 
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(HU) (see jChaSco :90j and references therein). Assume for simplicity that the Hamilto- 
nian is autonomous, that is, H = H(q,p). Now, first expands the generating function 
S2(t,q ,pi) as: 

oo 

^(t, go, pi) = gopi + ^i l Gi(go,pi), 

inserts expression into Hamilton-Jacobi equation Q2U|) and compares equal powers of 
t. This yields 

Gl(%,Pi) = H(q ,pi) 

1 (dH dH 

&2(go,Pi) : 



2 \dq£ dpiA 



do 

G 1 / d 2 H dH dH + d 2 H dH dH + d 2 H dH dH 

" "' ' 6 \dp 1A dpiB dq£ dq^ dpi A dq^ dq£ dp\ B dq^dq^ dp iA dpi B 



Using the truncated series, we obtain an approximated generating function: 

r 

S$(qk,Pk+i) = Qk ■ Pk+i + ^2 h r Gi(qk,Pk+i) 

i=i 

which defines a symplectic method of order r. 

Other approaches are also admissible without using higher derivatives of the Hamil- 
tonian H, for instance, symplectic or symplectic partitioned Runge-Kutta methods (see 
|HahiiWa:02[ISariCal : 94p. 



5 Discrete Hamiltonian systems 

In |ErbYan:92| Erbe and Yan have considered discrete linear Hamiltonian systems of 
the form: 

Ay(t) = B(t)y(t + l) + C(t)z(t) 
Az(t) = -A(t)y{t + l) - B T {t)z(t) 



where A, C are symmetric and / — B is invertible. Here Ay(t) = y(t + 1) — y(t), 
Az(t) = z(t + 1) - z(t) and y, z G R d . 

This problem is a particular case of a discrete Hamiltonian systems of the form 

Ay(t) = H z (t,y(t + l),z(t)) (33) 
Az(t) = -H y (t,y(t + l),z(t)) (34) 

where H(t, y, z) = ^(y T , z T ) ( J^*^ ^(t) ) ( z ) ' The s y m P lecticit y of the dis- 
crete linear Hamiltonian system was fully studied (see |ErbYan:92j . for instance, and 
references therein). The existence of a corresponding symplectic structure for discrete 
nonlinear Hamiltonian systems given by ()33|) and (|34j) was proposed by Ahlbrandt as 
an open problem ( Ahlb:9lH] and also Shi:02 ). 
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From the point of view of section |SJ this open problem is easily solved considering 
as generating function of the second kind the following one: 

St^ivit + 1), *(*)) = <t)y{t + 1) - H(t, y(t + 1), z(t)) . 
Then Eqs. (|33|) and ffity are precisely 

f y(t) = 9 -^(y(t + l),z(t)) 
\ z{t + l) = fi8^>(y(t + !),*(*)), 

which guarantees the symplecticity of the discrete Hamiltonian system. In order to 
find the canonical transformation associated to this generating map it is only necessary 
to impose the local condition (see |Arn:78j ): 

det ( d^st 1 ' t) (y(t + i),z(t) )_ LQ 

I dydz J 

Then, in a neighbourhood of a point satisfying the above condition, there exists a 
symplectomorphism defined by Eqs. (|3*3*|) and (|34|). 



Acknowledgments 

This work has been supported by grant BFM2001-2272 (Ministry of Science and Tech- 
nology, Spain). A. Santamarfa Merino wishes to thank the Programa de Formacion 
de Investigadores of the Departamento de Education, Universidades e Investigation of 
the Basque Government (Spain) for financial support. 



References 

[Ahlb:93] Ahlbrandt C D 1993 Equivalence of Discrete Euler Equations and Discrete 
Hamiltonian Systems, J. Math. Anal. Appl. 180, 498-478 

[Arn:78] Arnold V I 1978 Mathematical Methods of Classical Mechanics (Graduate 
Text in Mathematics 60, Springer- Verlag New York) 

[BobSus:99a] Bobenko A I and Y B Suris 1999 Discrete Lagrangian reduction, discrete 
Euler-Poincare equations, and semidirect products Lett. Math. Phys. 49, 79-93 

[BobSus:99b] Bobenko A I and Suris Y B 1999 Discrete time Lagrangian mechanics 
on Lie groups, with an application to the Lagrange top Comm. Math. Phys. 204 
147-188 

[Cadz:70] Cadzow J A 1970 Discrete calculus of variations Intern. J. Control. 11, 
393-407 

[ChaSco:90] Channell P J and Scovel C 1990 Symplectic integration of Hamiltonian 
Systems, Nonlinearity 3, 231-259 

[ChLeMa:94] Chinea D, de Leon M and Marrero J C 1994 The constraint algorithm 
for time-dependent Lagrangians, J. Math. Phys. 35 (7), 3410-3447 

[Cort:02] Cortes J 2002 Geometric, control and numerical aspects of nonholonomic 
systems (Lecture Notes in Mathematics, vol. 1793, Springer- Verlag) 



15 



[CorMar:01] Cortes J and Martinez S 2001 Nonholonomic integrators Nonlinearity 14, 
1365-1392 

[Dirac:64] Dirac P A M 1964 Lecture on Quantum Mechanics (Belfer Graduate School 
of Science, Yeshiva University, New York) 

[ErbYan:92] Erbe L H and Yan P 1992 Disconjugancy for linear Hamiltonian difference 
systems, J. Math. Anal. Appl. 167, 355-367 

[GotNes:79] Gotay M J and Nester J M 1979 Presymplectic Lagrangian systems I: The 
Constraint Algorithm and the Equivalence Theorem, Ann. Inst. Henri Poincare, 
A30 129-142 

[HaLuWa:02] Hairer E Lubich C and Wanner G 2002 Geometric Numerical In- 
tegration, Structure-Preserving Algorithms for Ordinary Differential Equations 
(Springer Series in Computational Mathematics 31, Springer- Verlag Berlin Hei- 
delberg) 

[JarNor:97a] Jaroszkiewicz G and Norton K 1997 Principes of discrete time mechanics 
I: Particle systems J. Phys. A 30, 3115-3144 

[JarNor:97b] Jaroszkiewicz G and Norton K 1997 Principes of discrete time mechanics 
II: Classical field theory J. Phys. A 30, 3145-3163 

[JorPol:64] Jordan B W and Polak E 1964 Theory of a class of discrete optimal control 
systems J. Electron. Control 17, 697-711 

[KaMaOr:99] Kane C Marsden J E and Ortiz M 1999 Symplectic energy-momentum 
integrators J. Math. Phys. 40, 3353-3371 

[KMOW:00] Kane C Marsden J E Ortiz M and West M 2000 Variational integrators 
and the Newmark algorithm for conservative and dissipative mechanical systems 
Internat. J. Numer. Math. Eng. 49, 1295-1325. 

[Lee:83] Lee T D 1983 Can time be a discrete dynamical variable? Phys. Lett., 122B, 
217-220 

[Lee:87] Lee T D 1987 Difference equations and conservation laws J. Statis. Phys., 46, 
843-860 

[LeMaMD:96] de Leon M Marrero J C and Martin de Diego D 1996 Time-dependent 
constrained Hamiltonian systems and Dirac brackets Journal of Physics A : Math. 
Gen. 29 6843-6859 

[LeMDSa:02a] de Leon M Martin de Diego D and Santamaria A 2002 Geometric in- 
tegrators and nonholonomic mechanics, Preprint IMAFF-CSIC 

[LeMDSa:02b] de Leon M, Martin de Diego D and Santamaria A 2003 Geometric 
numerical integration of nonholonomic systems and optimal control problems 2nd 
IFA C Workshop on lagrangian and Hamiltonian Methods for Nonlinear Control, 
Seville 2003, 163-168. 

[LeoMdD:2002] de Leon M and Martin de Diego D 2002 Variational integrators and 
time-dependent Lagrangian systems Rep. on Math. Phys 49 2/3, 183-192 

[Lew:86] Lewis F.L 1986 Optimal Control (John Wiley& Sons, New York) 

[Mae:80] Maeda S 1980 Canonical structure and symmetries for discrete systems Math. 
Japonica 25, 405-420 



16 



[Mae:81] Maeda S 1981 Extension of discrete Noether theorem Math. Japonica 26, 
85-90 

[MaPaSh:98] Marsden J E Patrick G W and Shkoller S 1998 Multisymplectic geometry, 
variational integrators, and nonlinear PDEs' Comm. Math. Phys. 199, 351- 395 

[MarWes:01] Marsden J E and West M 2001 Discrete mechanics and variational inte- 
grators Acta Numerica , 357-514 

[MosVes:91] Moser J and Veselov A P 1991 Discrete versions of some classical inte- 
grable systems and factorization of matrix polynomials Comm. Math. Phys. 139, 
217-243 

[Nijvan:90] Nijmeijer H and van der Schaft A J 1990 Nonlinear dynamical control 
systems (Springer- Verlag, New York) 

[NorJar:98] Norton K and Jaroszkiewicz G 1998 Principes of discrete time mechanics, 
III: Quantum field theory J. Phys. A 31, 977-1000 

[PKMO:02] Pandolfi A Kane C Marsden J E and Ortiz M 2002 Time-discretized vari- 
ational formulation of nonsmooth frictional contact Int. J. Num. Methods in En- 
gineering 53, 1801-1829 

[SanCal:94] Sanz-Serna J M and Calvo M P 1994 Numerical Hamiltonian Problems 
(Chapmanfe Hall, London) 

[Shi:02] Shi Y 2002 Symplectic structure of Discrete Hamiltonian Systems, J. Math. 
Anal. Appl. 266, 472-478 



17 



